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Abstract. This paper studies finite volume schemes for scalar hyperboHc conservation laws 
on evolving hypersurfaces of R^. We compare theoretical schemes assuming knowledge of 
all geometric quantities to (practical) schemes defined on moving polyhedra approximating 
the surface. For the former schemes error estimates have already been proven, but the 
implementation of such schemes is not feasible for complex geometries. The latter schemes, 
in contrast, only require (easily) computable geometric quantities and are thus more useful 
for actual computations. We prove that the difference between approximate solutions defined 
by the respective families of schemes is of the order of the mesh width. In particular, the 
practical scheme converges to the entropy solution with the same rate as the theoretical one. 
Numerical experiments show that the proven order of convergence is optimal. 



1. Introduction 

Hyperbolic conservation laws serve as models for a wide variety of applications in contin- 
uum dynamics. In many applications the physical domains of these problems are stationary 
or moving hypersurfaces. Examples of the former are in particular geophysical problems |27] 
and magnetohydrodynamics in the tachocline of the sun |15l I24j . Examples of the latter in- 
clude transport processes on cell surfaces [22] , surfactant flow on interfaces in multiphase flow 
[5] and petrol flow on a time dependent water surface. There are several recent approaches 
to the numerical computation of such equations. Numerical schemes for the shallow water 
equations on a rotating sphere can be found in [6l [TBI 123) . For the simulation of surfactant 
flow on interfaces we refer to [H HI [T7] . As we are interested in numerical analysis we focus 
on nonlinear scalar conservation laws as a model for these systems. The intense study of 
conservation laws posed on fixed Riemannian manifolds started within the last years. There 
are results about well-posedness [SJ [TUl [20] of the differential equations and about the conver- 
gence of appropriate finite volume schemes [21 [T3l [T^ [T9] . For recent developments on finite 
volume schemes for parabolic equations we refer to |21j . 
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In the previous error analysis for finite volume schemes approximating nonlinear conser- 
vation laws on manifolds the schemes were defined on curved elements lying on the curved 
surface and it was assumed that geometric quantities like lengths, areas and conormals are 
known exactly. While this is a reasonable assumption for schemes defined on general Rie- 
mannian manifolds or even more general structures ^LS^ with no ambient space, most engi- 
neering applications involve equations on hypersurfaces of and one aims at computing the 
geometry with the least effort. This is in particular important for moving surfaces where 
the geometric quantities have to be computed in each time step. Now the question arises to 
which extent an approximation of the geometry influences the order of convergence of the 
scheme. 

We consider the following initial value problem, posed on a family of closed, smooth 
hypersurfaces F = T{t) C M'^. For a derivation cf. [Ill I25j . For some T > 0, find 
u:Gt:= Ute[o,T] ^(i) x {t} ^ M with 

(1.1) ii + uVr-v + Vr-iiu,-,-) = in Gt, 

(1.2) n(-,0)=uo onF(O), 

where v is the velocity of the material points of the surface and uq : r(0) — )• M are initial 
data. For every u G M, t £ [0, T] the flux ¥{u, •, t) is a smooth vector field tangential to F(t), 
which depends Lipschitz on u and smoothly on t. Moreover, we impose the following growth 
condition 

(1.3) \Vr ■ f{u,x,t)\ <c + c\u\ yueR,{x,t) eGr 

for some constant c > 0. By u we denote the material derivative of u which is given by 

u{^t{x),t) := j^u{^t{x),t), 

where $t : F(0) — )• F(t) is a family of diffeomorphisms depending smoothly on t, such that ^"0 
is the identity on F(0). Obviously this excludes changes of the topology of F. We will assume 
that the movement of the surface and also the family $f is prescribed. A main result of this 
paper is a bound for the difference between two approximations of u. In particular, we will 
give an estimate for the difference between the flat approximate and the curved approximate 
solution. By curved approximate solution we refer to a numerical solution given by a finite 



volume scheme defined on the curved surface, cf. Section 2.2 , and by flat approximate solution 



we refer to a numerical solution given by a finite volume scheme defined on a polyhedron 



approximating the surface, cf. Section 2.3 We will see that the arising geometry errors can 
be neglected compared to the error between the curved approximate solution and the exact 
solution, i.e. both approximate solutions converge to the entropy solution with the same 
convergence rate. We will present numerical examples showing that the proven convergence 
rate is optimal under the assumptions for the numerical analysis. However, for most numerical 
experiments we observe higher orders of convergence. 

Our analysis also indicates that the geometry error poses an obstacle to the construction 
of higher order schemes. To this end we perform numerical experiments underlining in which 
manner the order of convergence of the higher order scheme is restricted by the approximation 
of the geometry. This shows that to obtain higher order convergence also the geometry of 
the manifold has to be approximated more accurately, cf. [9] in a finite element context. 



Geometric Error of FV Schemes for Conservation Laws on Evolving Surfaces 



3 



The outline of this paper is as foUows. In Section [2] we review the definition of finite volume 
schemes on moving curved surfaces and define finite volume schemes on moving polyhedra 
approximating the surfaces. The approximation errors for geometric quantities are established 
in Section [3j Section [4] is devoted to estimating the difference between the curved and the 
flat approximate solution. Finally, numerical experiments are given in Section [5| 

2. The Finite Volume Schemes 

This section is devoted to the construction of a family of triangulations Th{t) of the surfaces 
suitably linked to polyhedral approximations Th{t) of the surfaces. Afterwards we will recall 
the definition of a finite volume scheme on Th{t) which was considered in the hitherto error 
analysis and define a finite volume scheme on Th{t) which is an algorithm only relying on 
easily computable quantities. We mention that our triangulation as well as the definition of 
the finite volume scheme on Th is in the same spirit as the one Lenz et al. [21] used for the 
diffusion equation on evolving surfaces. 

2.1. Triangulation. We start by mentioning that there are neighbourhoods M{t) C M'^ of 
r(t) such that for every x S M{t) there is a unique point o(x, t) £ T{t) such that 

(2.1) X = a{x, t) + d{x, t)i/r(t)(a(x, t)), 

where d{-,t) denotes the signed distance function to r{t) and ^'r(t)(a(x, t)) the unit normal 
vector to r(t) pointing towards the non-compact component of M^\r(t). See jl2] for example. 

Let us choose a polyhedral surface Th{0) C A^(0) which consists of flat triangles such that 
the vertices of r/i(0) lie on r(0), and h is the length of the longest edge of Th{0). In addition 
we impose that the restriction of a|rh(o) • ^h{^) — ^ r(0) is one-to-one. We define Th{t) as the 
polyhedral surface that is constructed by moving the vertices of (t) via the diffeomorphism 

and connecting them with straight lines such that all triangulations share the same grid 
topology. A triangulation Th{t) of Th{t) is automatically given by the decomposition into 
faces. We define the triangulation Th{t) on T{t) as the image of Thit) under a(-, t)|r^^((). 
We will denote the curved cells with K{t) and the curved faces with e{t). A fiat quantity 
corresponding to some curved quantity is denoted by the same letter and a bar, e.g. let 
e{t) C r(t) be a curved face then e{t) = (a(-, t)|r^(j))~"^(e(t)). In order to refiect the fact that 
all triangulations share the same grid topology we introduce the following misuse of notation. 
We denote by K the family of all curved triangles relating to the same triangle -^(0) on 
Th{0). We do the same for e,K,e. Analogously by Th we denote the family of such families 
of triangles K. 

For later use we state the following Lemma summarizing geometric properties, whose 
derivation can be found in [T2] . 

Lemma 2.1. Let Th{t) be a polyhedral approximation ofT{t) as described above then there 
exists C = C{T) such that for all t € [0, T] 

(1) vr{t) = ^d{-,t), 

(2) ||d(.,t)|r,(t)IU-(r,{t)) <C/i2 . 

We will use the following notation. By hx_i^t) •= diam(iir(t)) we denote the diameter of 
each cell, furthermore h := max^g^.T] ™ax^(t) ^_ft:(t) and , \dK{t)\ are the Hausdorff 
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measures of K{t) and the boundary of K{t) respectively. When we write e{t) C dK{t) we 
mean e{t) to be a face of K{t). 

We need to impose the following assumption uniformly on all triangulations 7h{t). There 
is a constant number a > such that for each flat cell K{t) £ Thit) we have 

ahj^U) < \K{t)\ , 
(2.2) ^^^> ' ' 

a\dK{t)\ < hfi^^y 



Later on, we will see that (2.2) implies the respective estimate for the curved triangulation, 
cf. Remark 3.4 A consequence of ([2^ is that 2a^hp^,^^^ is a lower bound of the radius of 
the inner circle of K{t), which implies that the sizes of the angles in K{t) are bounded from 
below. Furthermore we denote by K{x,t) the supremum of the spectral norm of V^'r(t)(2^)- 
By straightforward continuity and compactness arguments n is uniformly bounded in space 
and time. 

2.2. The Finite Volume Scheme on Curved Elements. In this section we will briefly 
review the notion of finite volume schemes on moving curved surfaces. We consider a sequence 
of times = to < < *2 < • • • and set In '■= [tn,tn+i]- Moreover we assign to each n G N 
and K £ Th the term approximating the mean value of u on Ute/„ I^i^) ^ {^} and to each 
K £ Th and face e C dK a numerical flux function ^ : — t- M, which should approximate 



(2.3) f f {f{u{x,t),x,t),iXKAx,t))de{t)dt, 

Jin Je{t) 

where de{t) is the line element, ^K,e{x, t) is the unit conormal to e{t) pointing outwards from 
K{t) and (•, •) is the standard Euclidean inner product. Please note that ^K,e{t) is tangential 
to V{t). Then the finite volume scheme is given by 

Uk ■■= / no(x)dr(0), 
Jk(q) 

(2.4) ,,ra+l ._ \K{tn)\ n l-^"! \„(f \\ fu ( n n \ 



u'^ix, t) := for t G [t„, tn+i),x G K{t), 



where denotes the cell sharing face e with K and (ir(0) is the surface element. For 
the convergence analysis it was usually assumed [131 HH] that the used numerical fluxes are 
uniformly Lipschitz, consistent, conservative and monotone. Additionally, the CFL condition 

(2.5) tn+i -tn< 

oh 

has to be imposed to ensure stability, where L is the Lipschitz constant of the numerical 
fluxes. Lax-Friedrichs fluxes satisfying this condition are usually defined by 

(2.6) f'l.Au^v) := -f- \i {f{u,;t) + f{v,;t),fiKAt))deit)dt + A(u - v), 



In 2 Je{t) 



where A = ^Pu/H oo is an artificial viscosity coefficient ensuring the monotonicity of ^ and 
stabilizing the scheme. 
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2.3. The Finite Volume Scheme on Flat Elements. In this section we define a finite 
volume scheme on Th which is in the same spirit as (2.4) but only relies on easily accessible 
geometrical information. We assume that / is smoothly extended from Gt to the whole 
of Ute[oT]-^(^) ^ {^}- want to point out that the calculation of areas and lengths is 
straightforward for flat elements. As well, the approximation of integrals can be achieved 
using quadrature formulas by mapping cells and edges to a standard triangle and the unit 
interval, respectively, using affine linear maps. In this fashion we obtain for every time 
t £ [0,T] quadrature operators Qf^^^t^ : C^{K{t)) M, and (5e(i) : C'°(e(*)) ^ K of order 
Pi,P2 £ N, respectively. In addition for any compact interval / C [0, T] the term Qj : C^{I) — >■ 
M denotes a quadrature operator of order £ N. 

Before we can use the quadrature operators to define numerical fluxes we need to determine 
the "discrete" conormals. To each flat triangle K{t) we fix a unit normal j^^(i) ^y imposing 

(2.7) {i?K(t),'^miy)) >o, 

where y is the barycentre of K{t). We will see in Lemma 3.2 that '^x(t) converges to '^r(t){y) 
for /i — 7- 0. To each face e{t) and adjacent cell K{t) there is a unique unit tangent vector 
'^K(t),e(t) such that i'K(^t) ^ '^K{t),e(t) & conormal to e{t) pointing outward from K{t). Hence 



However in general 



this vector product is one candidate for ^J'K(t)^e{t)- 
(2.8) 

such that a choice like 

f^K(t),e(t) 

would lead to a loss of conservativity of the resulting numerical fluxes. Therefore we choose 
We define a numerical Lax-Friedrichs fiux and a finite- volume scheme: 



^K{t) X ^K{t),e{t) 



1 



1 



2\e( 



[{f{u,;-) + f{v, 



./^i?(.),e-(.)/ 

+ \iu 



(2.9) 



-,n+l 



\K{tn)\ 



\In\ 



K 



\K{tn+l)\ ^ |i^(t„+l 

u^(x,t) ■.= u\, for f G [t„,t„+i). 



' e^dK 
X G Kit), 



for some sufficiently large A > 0. Note that by (2.9)4 the function u is defined on Gt- 



3. Geometrical Estimates 

In this section we derive estimates for the approximation errors of the geometric quantities. 
Throughout this section we suppress the time dependence of all quantities. All the estimates 
can be derived uniformly in time. To obtain the geometrical estimates, we introduce the 
following lift operator. 
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Definition 3.1. Let U C Th. and g a function on U then we define a function on a\r,^{U) 
as 

Similarly we define the inverse of this lift operator by 

.-I 
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for a function g defined on some U gT. 



We begin our investigation with the differences between the normal vectors of the flat and 
curved elements. 

Lemma 3.2. There is a constant C such that for all flat cells K and every y £ K we have 



(3.1) 



< Ch. 



The constant C depends on derivatives of d, in particular on k. 

Proof. WLOG we can assume that ^ is a subset of {{x,y,0) E | y < 0} such that e = 
{(s, 0, 0) G I s G [0, he]} is one of its faces and (^'r)3 '(y) > for some y £ K. We start by 
showing that there exists some constant C > such that 

(3.2) \{i'r)i\<Ch, fori = 1,2. 

We recall that z^r = Vd, where d is the signed distance function to F. As the vertices of Th 
lie on r we know that there exists {x, y,0) £ K such that 

d{0, 0, 0) = 0, d{he, 0, 0) = 0, d{x, y, 0) = 0. 

Hence, the directional derivatives of d with respect to {x,y,0) and (1,0,0) need to van- 
ish somewhere in K. Thus their absolute value is of order 0{h) on K. Due to the angle 
condition ( |2.2[ ) an analogous inequality also holds for the directional derivative of d with 
respect to (0,1,0). As the directional derivative of d with respect to (1,0,0), (0,1,0) co- 
incides with (i^r)i, (1^)2 ) respectively, this proves (3.2). This immediately implies (i^r)3 = 
zby/l - Oih"^) = ±1 + 0{h'^). By assumption {1^)3 



have i/j^ = (0,0, 1) which proves (3.1). 



1 + 0{h ) everywhere and by (2.7) we 

□ 



Lemma 3.3. For the difference between the length of a curved edge e and the corresponding 
flat edge e we have 



(3.3) 



1 



< Ch^, 



and for the difference between the area of a curved cell K and the corresponding flat cell K 
we have 



(3.4) 



\K\ 
W\ 



1 



< Ch^, 



where C does not depend on h but on k. 

Furthermore let Cg be the parametrization of e over e given by a\e then we have 

(3.5) |||c;(s)|| - l| < C/i^. 
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Proof. We assume without loss of generality that K cM? x {0}. For small enough h we can 
parametrize the curved cell K according to (2.1 ) by a parametrization c = a\j^ : K ^ K CM.^ 
with 

c{xi,X2) = (xi,X2,0) - d{xi, X2,0)i^r{c{xi, X2)), 

where we suppressed the third coordinate in K. The ratio of volume elements of K and K 
with respect to the parametrization c is given by 



VWl '■= \/det(5r), 

where the matrix g is defined by 

For the parametrization c of -fC we have 

diC = ej — (Vd, Si) vr o c — d diC (Vtr)"^ o c for i = 1, 2, 

where Cj denotes the i-th standard unit vector. Due to the bounded curvature of F and 
Lemma I2. II we can show that 

(3.6) dic = ei- {ii^r)ii'r)oc + 0{h'^) fori = 1,2. 



Applying (3.1) we see that 

ur = ±(0, 0, 1) + (D{h) and {a, 1^) = (tr), = 0{h) for i = l,2. 
Thus, for the matrix g we have 

_ fl + 0{h^) 0{h^) 

l + 0(/l2) 

which implies for the volume element 

(3.7) dK = y^\dK = ^/lTO{h?)dK = dK + 0{h^)dK. 

Therefore, we arrive at 

^\-ldK < C\K\ 

K 



\K\ - \K\ 



for the error of the cell area which proves (3.4). 

To prove (Q and ([3^ we consider WLOG an edge e = {(s, 0, 0)|0 < s < /ig} C dK, 
where /ig denotes the length of e. The corresponding curved edge e is parametrized by 

(3.8) Ce{s) = c(s, 0) = (s, 0, 0) - d{s, 0, 0)z.r(ce(s)). 
Due to the bounded curvature of F we get for the derivative 

(3.9) 4(S) = (1,0,0) - UT{Ce{s)){uv)l{Ce{s))+0{h''). 



Applying (3.1) we get 
(3.10) 

and therefore 



\c',{s)\\ = l + 0{h' 



\c'{s) 



1 ds 



<C\e\ h\ 



□ 



Jan Giesselmann, Thomas Miiller 



Remark 3.4. Let us note that an analogous estimate to (2.2) for curved elements is an easy 
consequence of (2.2), (3.3), (3.4) and the fact \hj^ — hxl < C/i^, which is a consequence of 
Lemma 13.31 



Lemma 3.5. There is a constant C (depending on k) such that for all flat cells K, all flat 
edges e C OK and every x G e we have 



(3.11) 
(3.12) 
(3.13) 



(A 



X)) 



{l^K,e^l^KA^)) - 1 



< Ch, 



where t denotes a unit tangent vector to e. We want to point out that this estimate is 
independent of the sign oft. 



Proof. It is sufficient to show versions of (3.11 ) - (3.13) where flx,e substituted by xiK,e- 
Then analogous results for I'j^^^ x g are immediate. Indeed, estimates (3.11) - (3.13) follow 
because flf^ g is the mean of the vectors x t^^ g and x g. Firstly, we address the 
proof of (3.11 ). Let the same assumptions as in the proof of Lemma 3.2 hold and in addition 
let e be given by {{x, 0, 0) G M'^ | x G [0, /ig]}- We obviously have 



(3.14) 



%xt;^„- = (0,1,0). 



Note that the assumptions of the proof of Lemma 3.3 are satisfied. Hence we can use (3.9) 
i.e. the parametrization of e given by c satisfies 

(3.15) c'is) = (1,0, 0) - Mc{s)){ut)i{c{s)) + 0{h^), 



and coincides with t(c(s)) up to standardization. Hence, in view of (3.5) we obtain 
(3.16) t-\x) = (1,0,0) - Mcis)){ur)iic{s)) + 0{h^) 

for some s S [0,/ig]. Combining ( 3.14[ ) and (3.16) we find using (3.2) 

i^K X tK.,t-'{x)) = \{ur)2{c{smur)Ms))\ + 0{h^) < Ch\ 



which is (3.11). Concerning (3.12), 



< Ch 



holds because of (3.14) and (3.2). Thus, it remains to show ( [3.13 ). By definition t~'(j;), z^p'(x), /u^'g(x) 
form an orthonormal basis of M'^ and the vector x g is of unit length. This means that 
for every j; in e there exist bi{x),b2{x),bs{x) G M satisfying 6f(x) + 62(0;) + ^3(3;) = 1 such 
that 



(3.17) 



bi{x)t ^{x) + b2{x)u^\x) + b3{x)nj^\{x). 



We know from (3.11) and (3.12) that |6i(j;)|, |62(a;)| < Ch for some C > 0, which implies 
using Taylor expansion 

(3.18) bsix) = ±Vl + Oih^) = ±l + 0{h^). 
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Note that it only remains to show that in (3.18) the + holds. As 63 depends continuously 
on X it is sufficient to find one (a;i,0, 0) G K such that bs{xi) = 1 + 0{h?'). To that end we 
consider some xi, yi > such that 

7:(-yi,0] — > K, si-^(xi,s,0) 

is a curve leaving K through e. By definition the curve 7 given by 

7(s) := 7(s) - i/f '(7(s))fi(7(s)) 

is a curve in K leaving through e. This means we have 

(3.19) 0< (7'(0),M,^,e(7(0))). 



Due to (3.17), (3.18) and the fact that ^K,e is of unit length we already know that 

(3.20) ^i^,e = ±(0,l,0) + O(/i). 
We are able to compute 

(3.21) 7'(s) = (0, 1, 0) - (0, 1, Q){Vu^\^{s))fd{-f{s)) - ur{l{s)) {Ml{s)), (0, 1, 0)) 



(0,l,0) + O(/i), 



because Vi^p' is bounded, Lemma 2.1 and ( |3.2[ ). Inserting (3.20) and (3.21) in (3.19) we find 
(3.22) 0<±l + C'(/i), 



where ± is the sign from (3.18). Obviously for h sufficiently small (3.22) only holds for "+", 
which finishes the proof. □ 



4. Estimating the Difference Between Both Schemes 

This section is devoted to establishing a bound for the difference between the curved and 
fiat approximate solutions. To start with we investigate the difference between the numerical 
fiuxes defined on the flat and the curved triangulation respectively. 



Lemma 4.1. Let K, he some compact subset ofM?. Provided the quadrature operators Qe(t) 
and Qj^ are of order at least 1, then there is a constant C depending only on Gt and /C such 
that for the Lax- Friedrichs fiuxes (2.6) and (2.9)i with the same diffusion rate A the following 
inequality holds 



fleiu, v) - fl Ju, v) < Ch^ V (n, v)^K,, K^Th, dK. 
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Proof. We start by observing that the diffusive terms drop out, such that 



(4.1) 



In Je{t) 



{f{u, X, t),HK{t),e{t) {x))de{t) dt 



In Je{t) 



1 



-Ml 



{f{v, X, t),iJ,K{t),e{t) {x))de{t) dt 

|3^<3g(-)[(/(^''T)>/^i?(.),e(-))] 



As u and v appear symmetrically we will omit all terms containing the latter in our subsequent 
add zero several times 

fK,eiu,v)-f^^iu,v) 



analysis. We now add zero several times in (4.1) and get 
(4.2) 



T1 + T2 + T3 + T4 + T5 dt 



In 



with 



Ti{t) 
T2{t) 

Ut) 



{f{u,x,t),iiK{t),e{t){x))de{t)-f (/ {u,x,t),ii^,^.,Jx))de{t), 

e{t) J ^it) 



if {u,x,t),fij^,^.,Jx))de{t)-f if (n,x,t),^^(4)^g(j))de(t) 

e(t) ' ''e(t) 



(/ (-", X, t), g(i))de(t) - f {f{u, X, t), ^^(j) g(i))(ie(t) 

e(t) ■'e(t) 



iQi 



(/(n,x,t),/i^(t)^g(i))de(t) 



1 



iQi 



{f{u,x,-), ^^Ki■), si-)) M-) 



(/(n,x,-),/^i?(.),e(.))^e(-) - j^^Qe(.) [fiur,-),IJ'K {■),€{■)) 



In the following we will estimate the summands one by one. First, by properties of the 



quadrature operators Qi,^, Qe(t) and the CFL condition (2.5) 



(4.3) 



In 



Ti{t) dt 



as the integrands are sufficiently smooth. In particular, we use the fact that the surface 
evolves smoothly. Addressing the estimates for Ti,T2,T3 we will omit the time dependency 
as all three estimates are uniform in time. To establish an estimate for Ti we recall that we 



can parametrize e over e such that for the parametrisation Cg inequality (3.5) holds. We have 
(4.4) |ri|= /(r'(n,x),^^;,(x))(||c;(.)||-l)de- 



< 



where ||/||oo denotes the supremum of f{u,x,t) for {x,t) G Gt and n G /C. Next we turn to 
T3. Its estimate is based on the assumption that we have extended f{u, •) to M smoothly and 
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on the second statement of Lemma |2.1[ This leads to 
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(4.5) IT3I < 

This leaves T2. It is clear that 



f \u,x) - f{u,x) ||/i^ g|| de < C/i^. 



(4.6) 



IT2I < max 



Furthermore we find, as / is tangential to F, 



/ \U,X) = fi{u,x)t \x)+ f2{u,x)^JL^^{x), 



where t is a unit tangent vector to e and fi{u,x), f2{u,x) G M. Due to Lemma 3.5 we have 

(4.7) {f-'{u,x),fxjl^{x)) = f2{u,x), 

(4.8) {f-'{u,x),fiK,)=fi{u,x)0{h^) + f2{u,x) + f2(,u,x)0{h^). 



Obviously it holds |/i(u,x)|, \f2{u,x)\ < \\f\\oo such that inserting (4.7),( |48| ) into (4.6) gives 
(4.9) \T2\ < Ch^. 



Now the statement of the Lemma follows from (4.2) together with (4.3), (4.4), (4.5) and 
( 1491 ). n 

Our next step is to establish stability estimates for the curved and flat approximate solu- 
tion. Due to the geometry change of the surface F which might act as a source term we need 
the following lemma. 

Lemma 4.2. For every finite sequence of positive numbers {bn}n=i,...,N we have 

N 

(4.10) " I 



N I ^ h \ / N \ 

n(l + M< I + Et^ <exp 

71=1 V n=l / \n=l / 



Proof. From Jensen's inequality we know 
(4.11) 



f;in(l + 6„)<iVln(f;^) 

n=l \n=l / 



Applying the exponential function to (4.11) gives the first inequality in (4.10). The second 
inequality in (4.10) follows from the fact that 

N 



1 + 



N 



< exp(c) ViV G N, c G 



□ 

Now we can show a stability estimate for the curved scheme, the proof of which is mostly 
standard. 

Lemma 4.3. Provided the initial data satisfy uq G L°°(F(0)), then the solution of the curved 
scheme fulfils 



(4.12) 



i^+'l < (l + c|/„|)max{K|,max{KJ}} + c|/„| Vi^ gT^, 

eCoK 
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for some constant c and therefore 

(4.13) ||n''(t)||Lcx, < (IIuoIIl- +cr)exp(cr) VO < t < T. 



Proof. Invoking the consistency of the numerical flux functions we can rewrite (2.4) 



U 



71+1 

K 



\K{tn 



with 



CK,, 



"+1^1 V eCdK eCdK 

-\In\i / Vv f{u'k:X,t)dT{t)dt 

Jl„ JK{t„) y 

\In\ |efa)|/^,eK.^L)-/^,eK^^x) 



Due to the monotonicity of the numerical fluxes and the CFL condition (2.5) we have 

CK,e > 0, ^ CK,e < 1- 
eCdK 



Combining the growth condition (1.3) and the fact that \K{tn)\/\K{tn+i)\ < 1 + c|/„| we get 
(4.12) for a new, possibly larger constant c. Iteration of (4.12) implies 



n-l 



n-1 



71-1 



(4.14) max < JJ(1 + c|/fc|) max |n5^| + ^c|4| JJ(l + c|/j|). 

k=o ^^^'^ k=0 j=k+l 



Invoking (4.10) we obtain from (4.14) 



n-l 



max \u^\ < exp(cT)||no||L°° + V] exp(cr) < (||uo||l°° + cT) exp(cT). 

k=o 

□ 

As a technical ingredient for the stability estimate of the flat scheme and the error estimate 
we need the following lemma whose proof is given in the appendix. 

Lemma 4.4. For times tn,tn+i and corresponding cells K{tn), K(tn-\^i), K{tn), K{tn+i) the 
following estimates hold 

\K{tn)\ \K{tn+l] 



(4.15) 
(4.16) 



\K{tn)\ \K{tn+l 
\K{tn)\ \K{tn+l)\ 



\Kitn+l)\ \K{t, 



< Ch\tn+1 — tn\ 

< Ch\tn+1 — tn\ 



The stability estimate for the flat scheme is a combination of the stability estimate of the 
curved scheme and the estimate for the difference of the fluxes. 

Lemma 4.5. Provided the initial data satisfy uq E L°°(r(0)), then the solution of the fiat 
scheme fulfils 



(4.17) 



-,n+l\ 
''K I 



< (1 + 2(c + 1)|I„|) max{|72||, max{|M^J}} + 2(c + 1)|I„| + d\In\h 

eCdK 



Geometric Error of FV Schemes for Conservation Laws on Evolving Surfaces 



13 



for all K £ Th and < tn+i < T. Here c can be chosen as the same constant as in Lemma 
\4-3\ and d > is another constant. Therefore, for h sufficiently small, 

-.hi 



(4.18) 

where h := 2(c + 1). 
Proof. We have 



L°° < (IkolU- +6r)exp(6r) + 1, VO < t < T, 



u 



n+l 
K 



(4.19) 



\K{tn)\ 
'\K{tn+l) 

\K{tn)\ 
'\K{tn+l) 

eddK 



\ I V n^l g^g^ 

/ \^{'tn)\ ,n (-n -n \ _l^(^n)l_ Jn /-n -n \ 



We observe that because of (|4.16|) 

\K{tn)\ 



^ \K{tn)\ \K{tn)\ \K{tn+l)\ 

\K{tn+l)\ \K{tn+l)\\K{tn+l)\ \K{tn)\ 

< (1 + C |/„|) • (1 + C \In\ /l) < 1 + (C + 1) \In\ , 



(4.20) 



where c is the same constant as in Lemma |4.3[ for h small enough. Moreover, provided 
maxKeTh I^JI < A + 1 := (||uo||l°° + bT) exp(6T) + 1 we have 



(4.21) 



\Kitn)\ / \eitn)\ ,n 



\K{tn 



+1) 



\K{tn) 



Htn)\ 7n 



\K{tn)\ 



fRA^R^^lJ <Ch 



because of (3.3), ( |3.4[ ), and Lemma 4.1, Here we have used that for < A + 1, the 

numerical fluxes f^^{u,v), f^-{u,v) are uniformly bounded. 

Provided maxxeTh I'^^JI — ^ + 1 ^^'i h, sufficiently small, we obtain (4.17) by the same 
argumentation as in the proof of Lemma 4.3 for some d > 0. As obviously ||n'^(0) ||loo < ^ + 1 
we have by induction 

n—l n—l n—1 

(4.22) max\u],\<]l{l + b\h\)max\u%\ + ^ib\h\+dh\h\) J](l + 6|/,|) 

k=0 k=0 j=k+l 

^ (||''^o||l°° + bT) exp(6T) +dT exp{bT)h, 

^ V ' 

<A 



where b = 2(c+ 1). Equation (4.22) shows that our induction hypothesis, max^^-gT^ < 
A + 1, also holds for the next time step provided h < ^^^^^^^^^^ and t„ < T. This implies that 
(4.17) and (4.22) in fact hold for all tn < T. Thus, provided h is small enough, the assertion 
of the lemma follows by induction. □ 
In addition we need the fact that the curved scheme satisfies a discrete L^-contraction 
property. 
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Lemma 4.6. For given data and v^, let and , be defined according to the 

curved finite volume scheme (2.4). Then 



■+1)\\Uk 



n+1 n- 



K 



K 



The proof is analogous to the proof of the discrete L^-contraction property for finite volume 



schemes in Euclidean space, cf. [7j. We recall that the Lax-Friedrichs flux (2.6) satisfies the 
classical conservation and monotonicity conditions. 

For the difference between the curved and flat approximate solutions we obtain the follow- 
ing estimate. 



Theorem 4.7. Provided the quadrature operators Q^{t) '^''^d Qj^ are of order at least 1 for 
all t, n and the quadrature operators Qk(o) initial data uq are such that 



(4.23) 



\u''{0)-u\0)\\L^rm<Ch 



for some constant C. Then, for fixed T > 0, the error between the solution of the curved 
finite volume scheme (2.4) and the solution of the flat finite volume scheme (2.9) satisfies 



u''{T)-u^{T) 



< Ch. 



for some constant C depending on T, Gt, f, uq, provided the same diffusion rate A is used in 
both schemes. 



Remark 4.8. The curved approximate solution converges to the entropy solution of ( 1.1 )-( 1.2 ) 



with a convergence rate of cf. [T4|. Hence, invoking Theorem 4.7 the same kind of 

error bound holds for the flat approximate solution. 



of Theorem \4- 7[ Let n G N be such that T G then we have 



|n'^(r)-n'^(r)||^i(r) = E|K(t„+i) 



K 



u 



n+1 „-;n+l 



K 



U 



K 



E 

K 



\K{t^)\u^K-\In\ Yl Ie(tn)|/^,e«,<. 
eCdK 



|K(Wl)||j^(t n;i^,. 

|i^(Wi)l 



\K{t 



n+1) 



eCdK 



< i?l + i?2 + R3 + -R4 + R5 ) 
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where 



R,:=Y,\\K{tn)\ul-\In\ Yl Ie(tn)l/X,e«,<j 



K eCdK 

eCdK 
' \K{tn+l)\ 



U 



K 



K 



eCdK 



R,:=Y 



K 



1 



\K{t 
\K{t 



^) \In\ Y fKAul,U%) 
eCdK 



Because we consider the Lax-Friedrichs flux which is conservative and monotone and due to 



the CFL-condition (2.5), term Ri can be estimated via the discrete L^-contraction property 

Ri<Y\^^^- 



(Lemma 4.6 ) 



K 



The term R2 can be estimated using ( |4.15 ), we get 

\Kitn)\ \K{tn+l)\ 



(4.24) 



R2<YK\\^(*n) 



K 



\K{tn)\ \K{tn+l)\ 



< C\In\h. 



<C\I„\h 



Applying Lemma 3.3 and assumption ( |2.2| ) together with Remark 3.4 we get 

(4.25) ^3,^4<J^|/n| Y Ch^\fK,eiul,u]^j\<C\In\h. 



Based on Lemma 3.3 assumption ( |2.2[ ), Remark |3.4| and Lemma 4.1 we have 
(4.26) R5<CY l^nl Y ''^^ 

K eCdK 

Combining these estimates we thus obtain by iteration 



u^{T)-u''{T) 



Li(r) 



n+lj 



K 



U 



n+1 



K 



U 



K 



<^|K(t„)| \u^^-ul\+C\In\ h 
K 

<Y\KmK-u\\ + cTh 



K 



<C{T + l)h, 
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where the last step follows with (4.23). □ 



5. Numerical Experiments 

Numerical investigations based on the finite volume schemes defined in Section [2] are pre- 
sented in this section. The upshot of our experiments is three-fold. Firstly, under the present 
assumptions the order of convergence stated in Theorem |4.7| is optimal. This is demonstrated 
by Test Problem [T} Secondly, all of our experiments which include a sufficiently large nu- 



merical viscosity, i.e. A G ©(1) in (2.6), lead to a considerably higher experimental order of 
convergence (EOC) between 1 and 2 for the L^-difference between the flat and the curved 
approximate solution. Thirdly, the application of a finite volume scheme of second order to 
Test Problem [T] demonstrates that orders of convergence higher than 1 are not to be expected 
in general, if the geometry is not approximated sufficiently well, see Test Problem [5] In the 
following we will present several test cases. Thereafter, we will mention some implementation 
aspects. 

5.1. Test Problems. All test cases except Test Problem [7] use the geometrical setting Gt = 
X [0,1], i.e. T{t) = §2 for all t G [0,T], and T = 1. This is due to the fact, that we 
are able to compute the exact curved quantities only in this or similarly simple settings. In 
addition, let us fix the vector fields V{x) = p|(x2, — xi, 0)"^ and W{x) = p||(— X3, 0, xi)'^ for 

X = {xi,X2, X3) G §2. 

Test Problem 1 (u-independent flux function). We choose f = V as the flux function. Since 
/ neither depends on t nor on u and is divergence-free on any initial datum uq : — s- M 



is a stationary solution of the corresponding initial value problem (1.1)-(1.2). For initial 
values identically to zero the curved scheme conserves this stationary solution. Thus, the 
error between the curved and the flat approximate solution is equal to the error between the 
fiat approximate solution and the exact solution. The results for this test case for A = are 
plotted in Table [T} Note that due to duf = the numerical flux functions are monotone. 
This experiment shows, that under the assumptions from our convergence analysis 0{h) is 
indeed the optimal order of convergence. 

However, if we modify the numerical diffusion by setting A = vr in the numerical flux 
functions we achieve EOCs between 1 and 2 as can be seen in Table [Tj as well. 



Test Problem 2 (Advection across the poles). Let the flux function / be defined by f{u, x) = 
uW{x) for X G S^. Initial values are given by uo{x) = l{a;i>o.i5}(^)- In order to get monotone 
numerical fiux functions we set A = ^||9u/||oo = ^r. For this test case we obtain EOCs of 
almost 2, cf. Table [2} 

Test Problem 3 (Burgers along the latitudes). We choose a flux function of Burgers-type 
/ = f{u,x) = l/2u'^V{x) for X G and initial values uo{x) = l{xi>o.i5}(3;)- In order to get 
monotone numerical flux functions we set A = ^||5m/||oo = tt and obtain EOCs of almost 2, 
cf. Table H 

Test Problem 4 (Fully two-dimensional problem). In this test problem we consider a flux 
function / such that the corresponding initial value problem is not equivalent to a family 
of one-dimensional problems. Note that the flux functions from the previous test problems 
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Test Problem 


1 


A = 


Test Problem 


1 


A = 1 


level 


L^-difference 


EOC 


X^-difference 


EOC 





0.758314 






0.0119577 






1 


0.437173 




0.794 


0.0050082 




1.256 


2 


0.231999 




0.914 


0.0020877 




1.262 


3 


0.119190 




0.960 


0.0008286 




1.333 


4 


0.060372 




0.981 


0.0003137 




1.401 


5 


0.030378 




0.990 


0.0001165 




1.429 


6 


0.015237 




0.995 


0.0000439 




1.409 



Table 1. L -difference and EOCs between curved approximate solution 
(T) and flat approximate solution v!^ (T) from Test Problem [l] for differ- 
ent values A of numerical diffusion. 



Test Problem 2 



L^-difference I E"UC 



Test Problem 3 



L^-difference I E"UC 



Test Problem 4 



L^-difference I E"DC 



level 



0.112518 
0.039167 
0.011223 
0.002984 
0.000772 
0.000197 
0.000053 



1.523 
1.803 
1.911 
1.951 
1.967 
1.891 



0.0370831 
0.0133379 
0.0040350 
0.0011216 
0.0002992 
0.0000778 
0.0000199 



.475 
.725 
.847 
.906 
.944 
.966 



0.115867 
0.035202 
0.009566 
0.002475 
0.000630 
0.000159 
0.000040 



719 
880 
950 
975 
985 
989 



Table 2. L^-difference and EOCs between curved approximate solution 
^{T) and flat approximate solution u^{T) from Test Problems [2| [s] and |4j 



u 



have been of one-dimensional nature. To this end we define f{x,u) = uV{x) + l/2v?'W{x) 
for X G with initial values uq{x) = l|^^>o 15}(^) ^"^^ observe EOCs of almost 2, cf. Table 

m. 



Test Problem 5 (2nd order scheme applied to Test Problem [T]). The 

motivation of this test problem is to show that in general even higher order schemes, which 
are based on the flat finite volume schemes, are not able to achieve higher order convergence 
rates for smooth data. To this end, we apply a second order finite volume scheme (which is 
validated in Test Problem |6]) to Test Problem [T| This scheme is based on the flat finite volume 
scheme of first order (cf. Subsection 2.3) with A = enhanced with a linear reconstruction 
and a second order Runge-Kutta method for time evolution. In Table |3] we observe EOCs of 
almost 1. Indeed, the application of a second order finite volume scheme to Test Problem [T] 
gives exactly the same convergence rates as a first order scheme since the linear reconstruction 
on each cell does not affect the numerical fiux functions as / is independent of u. We like to 
point out that we do not have to compute the curved approximate solution as it coincides 
with the (constant) exact solution. 
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Test Problem 


5 




Test Problem 


6 




level 


L 


-difference 


EOC 


1/ -error 


EOC 







0.777427 








0.00362492 








1 




0.444068 


0.808 


0.00243102 


0.576 


2 




0.233521 


0.927 


0.00112593 


1.115 


3 




0.119553 


0.966 


0.00034768 


1.700 


4 




0.060461 


0.983 


0.00012006 


1.534 


5 




0.030400 


0.992 


0.00003713 


1.693 


6 




0.015242 


0.996 


0.00001116 


1.734 



Table 3. L -difference and EOCs between second order curved approximate 
solution (which equals the exact solution in this case) and a second order flat 
approximate solution from Test Problem [s] and L^-error between the exact 
solution from Test Problem [6] and its approximation by the second order finite 
volume scheme. 



Test Problem 6 (Validation of the 2nd order scheme) . This test problem serves as validation 
of the second order finite volume scheme. We consider smooth initial values 

mix) := ^ l{rix)<i}{x) exp 

with r{x) := ^^qj^^ and xq := (1,0, 0)^ and a flux function /(x, u) := uV{x), which transports 
the initial values around the sphere. For the error between the flat second order finite volume 
scheme (see Test Problem [5]) and the exact solution EOCs significantly higher than 1 are 
shown in Table [3l 

Test Problem 7 (Deforming Torus) . We consider a deforming torus as computational domain 
r and T = 4 as final time. Within the time interval [0, 2] the right half of the torus undergoes 
compression whereas the left half is stretched, while T(t) remains constant for t G [2,4]. 
We choose a Burgers-type flux function / = f{u,x) = ^u'^{x2, —xi,0)^ and constant initial 
values uq = 1. The time step size is chosen dynamically for each time step such that stability 
is guaranteed. In Figure [l] the numerical solution is shown at four different times. Note that 
in spite of the constant initial values, a shock wave is induced due to the change of geometry 
(compression and rarefaction) and the nonlinearity of the flux function. 

5.2. Implementation Aspects. 

5.2.1. Software. All simulations have been performed within the DUNE-FEM module which 
is based on the Distributed and Unified Numerics Environment (DUNE) |8]. As coarsest grid 
approximating the sphere we use an unstructured grid consisting of 632 triangles, see Figure 
[2] For finer computations we refine the coarse macro grid (level 0) and obtain up to 2.5 
million, triangles for the finest grid (level 6) whose vertices are projected onto the sphere, cf. 
Table H 



-2 (1 + r2(x)) 



(1 



{x)y 
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(c) t = 2.36. (d) t = 4. 

Figure 1. Flat approximate solution for Test Problem [t] for four different 
times. The computation was performed on a deforming polyhedron consisting 
of about 3 million, triangles. 

5.2.2. Exact Computation of Spherical Volume. For the curved finite volume scheme on the 
sphere the exact outer conormals, exact lengths of boundary segments and exact volumes of 
spherical triangles need to be computed. While the computation of the former two quantities 
is an easy geometric exercise, we use the formula from |26j for the computation of the latter. 

5.2.3. Exact Computation of Numerical Flux Functions. For the exact evaluation of the nu- 
merical flux function corresponding to an edge e of a grid cell K, quantities of the form 
j'^{V, fj,K,e) de have to be computed. Note that V can be written as V = x Vhy with 
hv{x) = 27rx3, where iy{x) := x denotes the outer unit normal to S^. As a result, similar to 
[To] , we deduce 



r {V, lXK,e) de= 4- {flK,e X i^, V/iy) de. 
e J e 
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level 


h 


size 





0.033664 


632 


1 


0.016833 


2528 


2 


0.008416 


10112 


3 


0.004208 


40448 


4 


0.002104 


161792 


5 


0.001052 


647168 


6 


0.000526 


2588672 



Table 4. Different refine- 
ment levels of the sphere 
grid. 




Figure 2. The sphere grid 
of level 0. 



As iJ,K,e X is a unit tangent vector to e, the integrand is a directional derivative along e and 
thus the integral can be computed by the evaluation at the endpoints of e. Obviously, the 
same applies to W with h]v{x) = 2t:x2 and W = u x Vhw- 

5.2.4. Computation of L^- Norms. We remark that the L^-differences between the flat end 
the curved approximate solutions are computed on the triangulation F/^. This does not have 
any influence on the convergence rates. 

Acknowledgements. We gratefully acknowledge that the work of Thomas Miiller was sup- 
ported by the German Research Foundation (DFG) via SFB TR 71 'Geometric Partial Dif- 
ferential Equations'. Jan Giesselmann would like to thank the German Research Foundation 
(DFG) for financial support of the project 'Modeling and sharp interface limits of local and 
non-local generalized Navier-Stokes-Korteweg Systems'. 



Appendix 



Here we give the proof of Lemma 4.4 
Proof. As we know that 

\Kitn)\ \K{tn)\ \K{tn)\ 

|^(Wi)l'|i^(Wi)l'|^(Wi)l " 
it is sufficient to prove 



\K{t, 



\Kitn 



(5.1) 



\Kit 



n+l) 



\K(t 



n+l J 



< C\In\h, 



2 + 0{h^), 



mtn)\' \K{tn)\ 

to obtain both assertions of the Lemma. 

Without loss of generality we assume the following situation: the triangle K{tn) is the 
convex hull of (0, 0, 0), {h, 0, 0), {x, y, 0). We define 

c^"(-, t) : r{tn) ^ m, ^^{., t) := ^■,t) o $(., tn)-\ 



Geometric Error of FV Schemes for Conservation Laws on Evolving Surfaces 21 

such that <^"(-, t„) is the identity mapping. In addition we define for vectors a, b G M'^ 



a0b G M^^'^aOb := 



(a, a) (a, b) 
(a,b) (b,b) 



and observe the identities 



(5.2) ||a X bf = det(a0 b) and det(a b) = det((a + Ab) b) VA G M. 

We denote the canonical projection K{tn) — ?■ K{tn) by c and abbreviate o c by 
Obviously we have 

(5.3) \Kitn)\ = \hy 
and 

4|if (t„+i)|2 =11 ($"((x, y, 0), t„+i) - cl>"((0, 0, 0), Wi)) 

X {^^{{h, 0, 0), Wi) - <I>"((0, 0, 0), Wi)) f 

= det ( ($"(c(x, y, 0), t„+i) - $"(c(0, 0, 0), t„+i)) 
($"(c(/i, 0, 0), t„+i) - $"(c(0, 0, 0), t„+i)) 



as the vertices of K{tn) and K{tn) coincide. Continuing the computation using (5.2) we get 



:det ^(x,y,0) + ^ (9,$"(c(x,2/,0),r)-a,$"(c(0,0,0),r)) dr 
(^h, 0, 0) + ^ (a,^"(c(/i, 0, 0), r) - 9,$"(c(0, 0, 0), r)) dr 
:det[ (^ix,y,0) + dydr'^'^iixuyi,0),T)y + d,dr<l>':iix2,y2,0),T)xdT 



(5.4) 

( (/i, 0, 0) + / (a,.a,$^((x3, 0, 0), r)/i) 



In 



det ^((0,y,0) + ^ aj,9,$^((xi,yi,0),r)ydr 
+ / {d^dr^^{{x2, y2, 0), r) - a,a,$^((x3, 0, 0), r)) X dr 
(^(/i,0,0) + ^ (5,5,$^((x3,0,0),r)/i)dr 
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Combining ib.M and (5.4) we get 



\Kit 



n+l) 



(5.5) 



:det (((0,1,0)+ / dydr<^>]^i{xi,yi,0),T)dT 
+ / (a,ar$^((x2,2/2,0),r)-5,5,$^((x3,0,0),r))-dr 

Jin y 



0( (1,0,0)+ / (a,5,$^((x3,0,0),T))(ir 

Jin 



=1 + (^(0, 1, 0), dydr<^':{{xuyi,0),T)dT 
+ /(0, 1, 0), / {d,dr^:{ix2,y2, 0), r) - 9,9,<J>^((x3, 0, 0), r)) ^ dr 

\ J In y 

+ (^(1, 0, 0), (a,a,$^((x3, 0, 0), r)) dT^ + 0(|I„|2). 



[,2x2 



Here we used that 

(5.6) det{A + B) = det A + trAiB) + 0{\\Bf), yA,B£ 
Now we turn to the quotient of the areas of the curved faces and remark that 

(5.7) |i^(Wl)l = \<^''{K{tn),tn+l)\+0{\In\h\K{tn+l)\), 

cf. f2T] for a proof. 



|<f"(i^(tn),t 



n+l) 



K(tn) 



(5.8) 



det ((9^^>J? (9y^>J? 
/det (a^$?0aj/$J 



det{dxC 5j,c) 



/?(tn) V det(9a;C0 5yc) 

{(x4,y4,0),tn+l)\K{tn)\. 



/det (9a; ^>J^ 



det(92^c dyc) 



Using (|3.6|) we get for the quotient 

in+l 



\K{tn)\^ det{dxcQdyC 

dxc + Jjdrdx^'^{{x4, y4, 0), r)dr) x(ayC + Jjdrdy'^'^iix^, y4, 0), T)dr 



(5.9) 



det {dxC 9yc) 
1 + (^(1,0,0),^ a,a,«>^((x4,y4,0),r)dr^ 

+ ^(0,1,0),^ a,9j,ci>^((x4,y4,0),T)dr^+O(|/„|/i). 
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The statement of the Lemma foUows by considering the difference of (5.9) and (5.5) and the 



fact that the smoothness of o c only depends on the smoothness of □ 
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